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Garnet-type  Li7La3Zr20i2  (LLZO)  is  a  solid  electrolyte  material  with  a  low-conductivity  tetrag¬ 
onal  and  a  high-conductivity  cubic  phase.  Using  density- functional  theory  and  variable  cell  shape 
molecular  dynamics  simulations,  we  show  that  the  tetragonal  phase  stability  is  dependent  on  a 
simultaneous  ordering  of  the  Li  ions  on  the  Li  sublattice  and  a  volume-preserving  tetragonal  dis¬ 
tortion  that  relieves  internal  structural  strain.  Supervalent  doping  introduces  vacancies  into  the  Li 
sublattice,  increasing  the  overall  entropy  and  reducing  the  free  energy  gain  from  ordering,  eventually 
stabilizing  the  cubic  phase.  We  show  that  the  critical  temperature  for  cubic  phase  stability  is  low¬ 
ered  as  Li  vacancy  concentration  (dopant  level)  is  raised  and  that  an  activated  hop  of  Li  ions  from 
one  crystallographic  site  to  another  always  accompanies  the  transition.  By  identifying  the  relevant 
mechanism  and  critical  concentrations  for  achieving  the  high  conductivity  phase,  this  work  shows 
how  targeted  synthesis  could  be  used  to  improve  electrolytic  performance. 


The  garnet  structure  material  Li7La3Zi'20i2  (LLZO) 
has  an  ionic  conductivity  that  varies  by  two  orders 
of  magnitude  depending  on  whether  synthesis  pro¬ 
duces  a  cubic  (crcubic=L9xl0-4  S/cm)  |j}  or  tetragonal 
(tftetra=l-63xl0~6  S/cm)  I  phase.  With  the  higher  con¬ 
ductivity,  LLZO  has  excellent  solid  electrolyte  material 
characteristics.  Its  stability  against  both  Li  metal  and 
standard  cathode  LiCo02  0,  combined  with  a  ~5  eV 
band  gap  [lj  and  high  ionic  conductivity,  make  it  suitable 
for  exploiting  the  full  voltage  difference  between  anode 
and  cathode  while  circumventing  the  safety  concerns  in¬ 
herent  to  all  current  liquid  electrolytes.  Initially,  the  pro¬ 
cess  by  which  the  cubic  ground  state  could  be  stabilized 
against  the  tetragonal  was  unknown  and  uncontrolled. 
The  relevant  factor  was  eventually  traced  down  to  the 
addition  of  A1  in  the  structure,  whether  via  accidental 
uptake  from  Al-containing  crucibles  Q  or  via  intentional 
doping  Si-  Inclusion  of  other  supervalent  ions,  in¬ 
cluding  Ta,  Nb,  and  Ga  j|  has  also  been  successful 
in  producing  the  high  conductivity  phase.  However,  the 
underlying  mechanism  that  controls  the  transition  has 
so  far  remained  a  mystery,  hindering  further  progress  to¬ 
ward  improving  the  material  for  practical  usage. 

Here  we  use  our  recently  developed  variable  cell  shape 
density-functional  theory  (DFT)  plus  molecular  dynam¬ 
ics  (MD)  method  to  investigate  the  driving  force  behind 
the  tetragonal  to  cubic  transition  which  subsequently 
raises  the  conductivity.  We  find  that  in  the  cubic  phase, 
the  Li  sublattice  is  always  disordered  (all  Li  symmetry 
sites  partially  occupied),  while  in  the  tetragonal  phase  it 
is  always  ordered  (all  Li  sites  either  full  or  empty).  We 
find  that  the  energy  is  lowered  by  a  simultaneous  order¬ 
ing  of  the  Li  atoms  that  relieves  Li— Li  Coulomb  repulsion 
but  unfavorably  distorts  the  ZrOe  octahedra  and  a  lattice 
distortion  that  restores  the  preferred  Zr— O  bonds.  The 
two  symmetry-breaking  but  volume-preserving  phenom¬ 


ena  always  occur  in  conjunction  and  either  one  alone  ac¬ 
tually  raises  the  total  energy  of  the  system.  When  Al3+  is 
doped  into  the  system,  charge  compensation  takes  place 
through  creation  of  Li+  vacancies  that  reduce  the  free  en¬ 
ergy  advantage  of  complete  ordering  on  the  Li  sublattice, 
eventually  leading  to  disorder  and  a  transition  to  cubic 
symmetry.  The  transition  is  always  signaled  by  a  sudden 
shift  of  Li  occupation  which  can  be  used  to  pinpoint  the 
critical  temperature  and  vacancy  concentration.  Using 
these  criteria,  we  estimate  that  the  critical  concentration 
of  A1  dopants  necessary  to  achieve  the  high-conductivity 
cubic  phase  of  Li7_2XAla;La3Zr20i2  is  x  =  0.2  (vacancy 
number  nva c  =  2x  =  0.4  per  formula  unit)  at  zero  tem¬ 
perature,  in  good  agreement  with  experiment.  We  show 
that  the  cubic  phase  can  be  reached  at  some  tempera¬ 
ture,  regardless  of  Li  content,  but  that  the  critical  tem¬ 
perature  drops  as  a  function  of  vacancy  number.  The 
understanding  uncovered  in  this  work  will  be  useful  for 
choosing  more  efficient  dopants  and  further  raising  the 
ionic  conductivity. 

The  Li  sublattices  in  the  cubic  and  tetragonal  LLZO 
phases  are  shown  in  Fig.  [T]  Li  positions  in  each  struc¬ 
ture  are  generally  referred  to  as  Li(l)  if  they  are  tetra- 
hedrally  coordinated  to  oxygen,  and  Li(2)  and  Li(3)  if 
they  are  octahedrally  coordinated.  To  avoid  confusion, 
we  refer  to  each  site  using  its  crystallographic  notation, 
with  superscript  c  or  t  to  designate  cubic  and  tetrag¬ 
onal  lattices,  respectively.  Each  tetrahedral  cubic  24dc 
site  is  surrounded  by  four  pairs  of  octahedrally  coordi¬ 
nated  96 hc  sites,  and  all  sites  are  partially  occupied  (see 
Fig.  [T|) .  Because  of  Coulomb  repulsion,  it  is  energetically 
prohibited  for  both  members  of  each  pair  of  96 hc  sites 
to  be  occupied,  and  if  a  particular  24dc  site  is  occupied, 
the  adjacent  96 hc  sites  are  consistently  unoccupied  Q. 
The  tetragonal  distortion  transforms  the  cubic  24dc  sites 
into  fully  occupied  8o*  sites,  often  denoted  as  Li(l),  and 
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FIG.  1:  (Color  online)  Li  sublattice  in  the  cubic  (left)  and 
tetragonal  (right)  phases  of  LLZO.  All  Li  positions  are  in¬ 
cluded,  although  in  the  cubic  phase  not  all  are  occupied.  The 
Li(l)  atoms  (80*  and  24dc)  are  large  gray  (gold  online),  Li(2) 
or  16/*  are  white,  and  Li(3)  or  32<?t  are  dark  gray.  The  cubic 
Li(l)  positions  that  become  vacant  upon  transition  to  the  or¬ 
dered  tetragonal  structure  (16e*)  are  indicated  by  small  (gold 
online)  spheres. 


unoccupied  16e*  sites.  The  96 hc  sites  are  transformed 
into  fully  occupied  16/*  and  32 <7*  sites,  often  denoted  as 
Li(2)  and  Li(3),  respectively,  with  the  rest  fully  unoccu¬ 
pied.  As  we  will  show,  a  shift  from  8a*  sites  to  16e*  sites, 
both  subsets  of  24dc,  always  accompanies  the  tetragonal 
to  cubic  transition. 


All  of  our  simulations  use  density-functional  theory, 
with  the  Perdew-Burke-Ernzerhoff  exchange-correlation 
functional  @,0  as  implemented  in  the  Vienna  Ab-initio 
Simulation  Package  (VASP)  projector-augmented- wave 
software  version  5.2.12  ll],[l2|.  The  calculations  used  an 
energy  cutoff  of  400  eV,  except  for  the  relaxations  used  to 
explain  the  microscopic  basis  of  the  structural  transition, 
which  used  a  cutoff  of  600  eV.  Molecular  dynamics  (MD) 
time  evolution  is  carried  out  using  the  velocity  Verlet  al¬ 
gorithm,  as  implemented  in  the  lib  Atoms  13]  package, 
with  a  time  step  of  1  fs  and  using  forces  and  stresses  cal¬ 
culated  at  each  time  step  with  VASP.  All  MD  simulations 
use  constant  temperature  and  stress.  The  time  propaga¬ 
tion  algorithm,  specified  in  Supplementary  Information, 
is  a  rediscretization,  using  the  ideas  from  Ref.  0  and 
Ref.  Il5,  of  the  Langevin  constant  pressure  algorithm  in 
Ref.  [16 .  for  a  modified  version  of  the  equations  of  motion 
of  Ref.  [lit  Each  simulation  starts  from  a  unit  cell  of  the 
experimental  tetragonal  structure  of  Ref.  0,  for  a  total  of 
8  formula  units.  The  structure  is  relaxed  with  respect  to 
ionic  positions  and  unit  cell  size  and  shape  using  VASP, 
and  then  simulated  at  finite  temperature  and  zero  stress. 
While  the  MD  simulation  results  implicitly  include  the 
effects  of  entropy,  all  energies  we  calculate  explicitly  and 
quote  here  include  only  DFT  total  internal  energy.  The 
system  is  allowed  to  evolve  for  at  least  24  ps  with  ionic 
motion  as  well  as  cell  shape  and  volume  changes,  so  it 
can  spontaneously  transform  into  other  structures.  To 
determine  the  effect  of  composition  and  temperature  we 


simulate  the  stoichiometric  system,  with  56  Li  atoms,  as 
well  as  systems  with  1,  2,  and  4  Li+  ion  vacancies  per 
simulated  cell,  at  temperatures  ranging  from  300  K  to 
1300  K.  The  net  negative  charge  of  the  vacancies  is  com¬ 
pensated  by  a  uniform  positive  background  charge.  The 
vacancies  are  formed  by  removing  Li"1"  ions  randomly  se¬ 
lected  from  the  tetragonal  16/*  and  32gl  sites,  where  the 
vacancy  formation  energy  is  about  100  meV  lower  than 
at  the  8a*  sites.  The  crystallographic  site  identity  of  each 
atom  is  determined  during  each  MD  trajectory  (see  Sup¬ 
plementary  Information). 

An  example  of  the  time  evolution  of  one  system,  with 
7XVac  =  0.25  and  T  =  600  K,  is  plotted  in  Fig.  [2]  The  ra¬ 
tios  of  the  lattice  constants  along  x  and  y  to  that  along  z 
initially  show  the  tetragonal  structure  with  one  axis  (az) 
about  3%  shorter  than  the  others.  The  system  trans¬ 
forms  into  a  cubic  phase  where  both  axis  ratios  fluctuate 
around  1  at  t  ~  5  ps.  It  fluctuates  back  to  a  tetragonal 
phase  at  t  ~  15  ps,  and  then  again  to  cubic  at  t  ~  30  ps 
where  it  remains  for  the  rest  of  the  simulation.  The  vol¬ 
ume  is  not  affected  by  the  unit  cell  shape  change.  We 
also  plot  the  occupancies  of  various  crystallographic  sites 
of  the  cubic  and  tetragonal  structures,  computed  as  de¬ 
scribed  above.  We  see  a  perfect  correlation  between  the 
symmetry  of  the  unit  cell  parameters  and  the  occupancy 
of  the  tetragonal  sites.  The  8a*  occupancy  is  initially 
near  1,  as  it  is  in  the  experimental  structure,  and  most 
of  the  remaining  atoms  are  identified  as  16/*  and  32 <?*. 
When  projected  into  the  cubic  structure  the  8a*  atoms 
are  identified  as  24<ic  and  the  16/*  and  32 gt  atoms  are 
identified  as  96hc,  as  expected  by  symmetry.  Whenever 
the  system  transforms  into  the  tetragonal  structure  the 
occupancies  of  8a*  and  16ft+32gt  sites  drop.  Since  we  do 
not  see  a  corresponding  change  in  24dc  and  96 hc,  we  con¬ 
clude  that  the  atoms  are  moving  from  the  subset  of  the 
cubic  sites  that  are  occupied  in  the  tetragonal  ordering 
to  the  other  cubic  structure  sites.  This  is  indeed  seen  in 
the  increased  occupancy  of  16e*  sites,  which  correspond 
to  24dc  sites  that  are  unoccupied  in  the  experimental 
tetragonal  structure. 

The  unit  cell  shape,  8a*  occupancy,  and  unit  cell  vol¬ 
ume,  averaged  over  the  last  5  ps  of  each  run,  are  listed 
in  Table  ffl  At  each  vacancy  concentration  the  system 
transforms  from  an  ordered  tetragonal  structure  at  low 
T  to  a  cubic  structure  at  higher  T,  and  the  transition 
temperature  goes  down  with  increasing  vacancy  concen¬ 
tration.  Note  that  at  nva c  =  0.5  and  T  =  300  K  our  sim¬ 
ulations  predict  an  ordered  tetragonal  structure  that  is 
different  from  the  experimentally  observed  low  T  tetrag¬ 
onal  structure  for  the  stoichiometric  composition.  The 
structure  we  find  has  one  long  and  two  short  axes,  an 
occupation  of  the  original  tetragonal  8a*  sites  of  0.5,  and 
an  occupation  of  1  of  the  original  tetragonal  16e*  sites 
(which  are  unoccupied  in  the  stoichiometric  tetragonal 
structure).  The  volume  is  much  more  sensitive  to  tem¬ 
perature  than  to  vacancy  number,  with  a  linear  thermal 
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FIG.  2:  (Color  online)  Evolution  over  time  of  structural  and 
site  occupation  quantities  for  a  sample  system  with  nv ac  =  2 
at  T  =  600  K.  Top:  unit  cell  shape  ( ax/az  blue,  ay/az  red) 
and  volume  (black).  Middle:  96hc  (black)  and  16/t+32grt 
(red)  lattice  site  occupation.  Bottom:  24dc  (black),  80*  (red) 
and  16e*  (blue  with  symbols)  lattice  site  occupations. 

expansion  coefficient  (from  a  linear  fit  of  F1/3  vs.  T  to 
the  nv ac  =0.5  results,  where  we  have  the  widest  range 
of  temperatures)  of  2.2  x  10-4  K-1.  The  vacancy  num¬ 
ber  dependence  (from  a  linear  fit  of  V  vs.  nvac  to  the 
T  =  600  K  and  T  =  800  K  results),  which  is  the  va¬ 
cancy  formation  volume,  is  10±1.5  A3 /vacancy.  Note 
that  this  value  is  dependent  on  the  charge  compensation 
mechanism,  which  is  a  uniform  background  charge  in  this 
calculation,  compared  to  any  of  a  variety  of  supervalent 
dopants  in  experiment.  In  all  cases,  a  transition  to  cu¬ 
bic  symmetry  is  preceded  by  a  sudden  decrease  in  the 
80*  occupation  and  any  return  to  tetragonal  symmetry 
is  characterized  by  a  reoccupation  of  these  sites. 

We  have  also  computed  the  energy  difference  between 
the  tetragonal  and  cubic  structures  as  a  function  of  va¬ 
cancy  concentration.  For  the  tetragonal  structure  we  re¬ 
laxed  10  structures  with  random  vacancies  at  16/*  and 
32 gt  sites.  For  the  cubic  structure,  which  is  disordered 
even  at  the  stoichiometric  composition,  we  use  10  con¬ 
figurations  from  a  uniform  sampling  of  all  configurations 
that  obey  the  restriction  that  no  pairs  of  nearest  neigh¬ 
bors  (adjacent  96 hc  pairs,  or  a  24dc  and  its  nearest  neigh¬ 
bor  96 hc)  are  simultaneously  occupied.  We  plot  the  en¬ 
ergy  difference  between  the  minima  and  means  of  each 
of  the  ten  structures  in  Fig.  [3]  The  tetragonal  structure 
is  energetically  favored  for  nvac  <  0.25,  and  the  cubic 
is  favored  for  ?zvac  >  0.5.  Entropy  effects  will  shift  the 
equivalent  curves  for  the  free  energies.  We  expect  that 
this  shift  will  reduce  the  free  energy  of  the  tetragonal 
phase  more  than  the  cubic  phase,  because  the  tetrago- 


FIG.  3:  (Color  online)  Energy  difference  between  tetragonal 
and  cubic  structures  as  a  function  of  vacancy  number,  for 
minimum  energy  configuration  (solid  red  line)  and  mean  con¬ 
figuration  energy  (dashed  blue  line). 
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FIG.  4:  Relaxed  energies  (closed  symbols)  of  10  disordered 
configurations  and  one  at  the  experimentally  observed  tetrag¬ 
onal  phase  order,  for  a  cubic  unit  cell  and  for  a  fully  relaxed 
one,  relative  to  the  lowest  energy  fully  relaxed  ordered  tetrag¬ 
onal  cell.  The  energy  of  a  model  including  only  Zr06  distor¬ 
tions  is  also  plotted  (open  symbol). 

nal  phase  starts  out  ordered  and  therefore  gains  more 
entropy  with  the  introduction  of  vacancies  than  the  al¬ 
ready  disordered  cubic  phase.  The  transition  vacancy 
concentration  will  therefore  shift  to  higher  values  as  the 
temperature  increases,  although  we  expect  this  shift  to 
be  small  at  T  =  300  K. 

To  understand  the  relationship  between  Li  or¬ 
der/disorder  and  the  tetragonal/cubic  lattice  parameters, 
we  perform  total  energy  calculations  for  the  stoichiomet¬ 
ric  system  in  10  disordered  configurations  (again  from  a 
uniform  sampling  of  configurations  obeying  first-neighbor 
exclusion)  and  in  the  experimentally  observed  Li  order  in 
the  tetragonal  phase.  We  first  relax  each  one  with  respect 
to  atomic  positions  and  unit  cell  volume  but  constrain  the 
unit  cell  shape  to  remain  cubic,  and  then  relax  with  full 
freedom  for  the  unit  cell  size  and  shape  as  well  as  atomic 
positions.  The  total  energies  are  plotted  in  Fig.  [I]  The 
energies  of  the  disordered  system  are  generally  highest, 
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TABLE  I:  Mean  unit  cell  distortion  ( aa/az  —  1  for  a=x  and  y,  in  %)  (top)  and  occupancy  of  8 a*  sites  per  formula  unit  (bottom), 
averaged  over  the  last  5  ps  of  each  run,  for  different  vacancy  numbers  nvac  and  temperatures.  Transition  temperature  Tc  is  also 
indicated  for  each  quantity  with  a  clear  signal  of  a  structural  transition  (see  text). 
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with  only  a  small  gain  (and  correspondingly  small  change 
in  unit  cell  shape)  from  relaxing  the  cubic  unit  cell  con¬ 
straint.  Imposing  the  Li  order  while  maintaining  a  cubic 
unit  cell  actually  increases  the  energy  slightly  as  com¬ 
pared  with  the  mean  disordered  energy.  Only  once  the 
ordered  cell  is  allowed  to  relax  to  the  tetragonal  shape 
does  its  energy  become  lower  than  the  disordered  system, 
and  the  ordered  system  emerges  as  the  T  =  OK  ground 
state. 

To  understand  the  source  of  the  coupling  between  en¬ 
ergy  and  structure,  we  calculated  the  pair  distribution 
functions  (PDF)  of  various  ion  pairs  in  the  system  for 
the  cubic  and  tetragonal  ordered  systems,  along  with  a 
simple  point  charge  model,  based  on  the  nominal  ionic 
charges.  The  point  charge  model  shows  that  the  overall 
Coulombic  energy  is  actually  increased  upon  relaxation 
to  a  tetragonal  unit  cell,  indicating  that  the  energy  gain 
is  elsewhere.  The  PDFs  show  that  the  La3+— 02~  dis¬ 
tances  do  not  vary,  presumably  since  these  are  the  two 
largest  ionic  charges  in  the  system.  Li— O  and  Li— Li  dis¬ 
tances  do  change,  but  the  highly  ionic  nature  of  these 
interactions  means  that  such  changes  are  nearly  entirely 
accounted  for  in  the  point  charge  model.  The  shifting  of 
Li— Li  and  Li— O  distances,  combined  with  the  rigidity 
of  the  La— O  spacing,  results  in  some  distortion  of  the 
ZrC>6  octahedra.  Unlike  the  other  pairs,  the  Zr— O  in¬ 
teraction  is  at  least  partially  covalent.  In  the  cubic  cell, 
the  Zr— O  bond  lengths  are  2.130±0.020  A  and  O— Zr— O 
bond  angles  are  180±4.0°.  The  relaxation  to  a  tetrag¬ 
onal  cell  relieves  this  distortion,  restoring  the  octahe¬ 
dra  to  a  uniform  Zr— O  bond  length  of  2.125±0.005  A 
and  a  O— Zr— O  bond  angle  of  180±0.01°.  To  estimate 
the  energetic  contribution  of  ZrC>6  octahedra  distortions 
we  parametrize  the  calculated  total  energies  of  the  or¬ 
dered  tetragonal  cell  as  a  quadratic  function  of  Zr— O 
bond  length  and  O— Zr— O  bond  angle  by  computing  the 
energies  for  small  displacements  of  an  O  atom.  This 


parametrization  predicts  an  energy  for  the  ordered  cubic 
structure  of  0.083  eV/formula  unit  relative  to  the  tetrag¬ 
onal  structure  (Fig.  H|),  very  nearly  equal  to  the  DFT 
energy  difference.  We  therefore  conclude  that  lithium 
ordering  to  relieve  Li-Li  Coulomb  repulsion  leads  to  in¬ 
ternal  rearrangements  of  ions  to  maintain  La3+— 02+  dis¬ 
tances,  which  leads  to  a  lattice  distortion  that  allows  the 
ZrC>6  octahedra  to  preserve  their  preferred  shape. 

In  summary,  our  variable  cell  shape  DFT  MD  sim¬ 
ulations  of  LLZO  show  that  the  DFT  ground  state  at 
low  temperatures  is  the  experimental  ordered  tetragonal 
structure,  and  that  at  higher  temperatures  the  system 
transforms  into  the  experimental  disordered  cubic  struc¬ 
ture.  The  transition  temperature  decreases  with  increas¬ 
ing  Li+  vacancy  concentration,  and  the  disordered  cubic 
structure  has  lower  energy  when  the  number  of  vacancies 
per  formula  unit  is  larger  than  about  0.4.  These  results 
are  in  agreement  with  the  experimental  observations  of  a 
transition  as  a  function  of  Li+  vacancy  concentration  . 
The  microscopic  cause  of  the  structural  transition  is  the 
coupling  between  the  unit  cell  shape  and  the  hopping  of 
atoms  from  the  subset  of  the  disordered  cubic  sites  occu¬ 
pied  in  the  ordered  tetragonal  structure  to  the  remain¬ 
ing  cubic  sites,  which  are  unoccupied  in  the  tetragonal 
structure.  The  relative  stability  of  the  ordered  tetrago¬ 
nal  low  temperature  structure  is  driven  by  the  ordering, 
which  reduces  Li— Li  Coulomb  repulsion  but  distorts  the 
ZrC>6  octahedra,  and  the  tetragonal  distortion  which  al¬ 
lows  these  octahedra  to  return  to  their  preferred  high- 
symmetry  shape.  Our  simulations  enable  us  to  explain 
the  atomistic  mechanism  behind  this  finite  temperature 
structural  transformation  in  a  complex  material,  and  pre¬ 
dict  the  number  of  vacancies  necessary  to  achieve  the 
high  conductivity  material.  This  should  enable  better 
doping  schemes  that  optimize  ionic  conductivity  by  pro¬ 
viding  the  requisite  number  of  vacancies  with  as  few  arti¬ 
ficial  dopants  as  possible,  thereby  realizing  the  potential 
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of  LLZO  as  a  truly  stable  solid  electrolyte  material. 
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